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We present the scalar leptoquark pair production process at the LHC computed at the next-to-leading order 
in QCD, matched to the Pythia parton shower using the MC@NLO formalism. We consider the leading order 
decay of a leptoquark to a lepton (e, p, r or and a jet and observe the effects of parton shower on 

the final states. For demonstration, we display the kinametical distributions of a selection of observables along 
with their scale uncertainties for the 13 TeV LHC. We also present pair production cross sections and if-factors 
using massless five-quark flavor scheme for different LHC center-of-mass energies. The complete stand-alone 
code is available online. 
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I. INTRODUCTION 


Leptoquark (LQ or iq) is the generic name of bosons 
(scalars or vectors) that can couple to a quark and a lepton 
simultaneously. They appear in different beyond the Stan¬ 
dard Model (BSM) scenarios like the models with quark 
lepton compositeness [1], Pati-Salam models [2], SU{5) 
grand unified theories [3], the colored Zee-Babu model [4], 
i?-parity violating supersymmetric models [5] etc. The 
LHC is looking for their signatures [6-16]. 

Recently, CMS has searched for the first two generations 
of scalar LQs at the 8 TeV LHC with 19.7 fb“^ of integrated 
luminosity in two different channels - ££jj and [15] 
(also see Refs. [12] and [13]) (we use £ to denote a charged 
lepton from the first two generations and £ for any charged 
lepton, i.e., £ = and £ = ,t^). Consider¬ 

ing only LQ pair production, the analysis puts 95% C.L. 
mass exclusion limits (ELs) on the first generation LQ at 
Mt^ = 1005 (845) GeV for /3 = 1 (0.5), where /3 is the 
branching fraction for a LQ to decay to an electron-quark 
pair and denotes the mass of LQs. For the second gen¬ 
eration, the corresponding limits are put at 1080 (760) GeV. 
CMS also excludes third generation LQs decaying to a r 
and a 6-quark with /3 = 1 up to 740 GeV [12]. In the first 
generation search, mild excesses of events compared to the 
Standard Model (SM) background in both channels for LQs 
with mass around 650 GeV were observed. Currently, these 
excesses have attracted considerable attention in the litera¬ 
ture [17-27]. 

For the first generation LQ, a more stringent limit comes 
from ATLAS. With 20 fb“^ of integrated luminosity at the 
8 TeV LHC, ATLAS rules out first (second) generation LQs 
decaying to an e (a p) and a jet with /3 = 1 up to 1050 
(1000) GeV [9]. For the third generation, ATLAS rules out 


LQs decaying to a 6-quark and a i/t with /3 = 0 up to 625 
GeV whereas for LQs decaying to a f-quark and a with 
/3 = 0 the excluded range is between 210 and 640 GeV [9]. 

Most of the experimental analyses use the fixed order 
(FO) result for the LQ pair production cross section com¬ 
puted at the next-to-leading order (NLO) of QCD [28]. At 
FO, 0{as) contributions consist of real QCD emissions 
from the Born subprocesses and the interference between 
the Born diagrams and the one-loop corrected Born dia¬ 
grams. However, to get a more realistic picture of the fi¬ 
nal state particles, one has to combine the FO NLO results 
with parton showers (PS) that consistently resum large log¬ 
arithmic terms in the collinear limit, thereby expanding the 
coverage of the kinematical region. In this paper, we com¬ 
pute the pair production of scalar LQs at the NLO of QCD 
accuracy; i.e., we analyze the 0{as) corrections to the fol¬ 
lowing leading order (LO) partonic subprocess at the LHC, 
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where the curved connections above or below mark a pair 
of a lepton {£ or ix = Vr) and a jet (j) coming from 

the decay of a LQ, retaining all the spin-correlation effects 
at the LO accuracy and then we match the fixed order re¬ 
sult with the Pythia PS. Computation of such multijet in¬ 
clusive event samples is important to properly interpret ex¬ 
perimental data.* Compared to FO calculations, matching 
with the PS improves the accuracy of theoretical estima¬ 
tions of various kinematic distributions. These improved 
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^ Not only the PS but inclusion of other processes with similar final states 
can also affect the ELs significantly — see Ref. [29] to see how inclu¬ 
sive single production events in the pair production signal can affect the 
ELs significantly. Combining single and pair productions in signal sim¬ 
ulations leads to a more realistic signal [30, 31] estimation in general. 
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FIG. 1: The LO lepton exchange diagram for pair production of 
first generation LQ. 


and more precise distributions can play a crucial role if ad¬ 
vanced techniques like multivariate analysis are employed 
to separate LQ signal from SM background or differenti¬ 
ate LQs from other BSM candidate with similar signatures 
(like leptogluons for example). 

While matching the FO NLO results with PS, we use the 
MC@NLO formalism [32] that avoids double counting in 
the hard and collinear regions at the NLO level by adding 
suitable counterterms. As our NLOh-PS computation is per¬ 
formed using the MadGraph5_aMC@NLO framework 
[33], it is flexible to any choice of cuts, scales, parton distri¬ 
bution functions (PDFs), etc. We use this paper to elaborate 
the details of the computation and demonstrate its reliabil¬ 
ity and applicability. Hence, it will facilitate the experi¬ 
mentalists to generate more realistic QCD NLOh-PS level 
events in future. 


II. NLO+PS COMPUTATION 


We use the same parametrization as in Ref. [29] and con¬ 
sider the following simplified models for scalar LQs with 
electromagnetic charges, —1/3 and 2/3 respectively. 


c/)K 
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where. 





A| = /3A2, Xl = {l-(3)X\ (4) 


assuming there is no unknown decay mode of the LQs and 
there is no mixing among different generations. Here pL 
and pft = (l~t?L) are the charged lepton chirality fractions, 
i.e., pl {vk) gives the fraction of charged leptons coming 
from a LQ decay that are left handed (right handed). We 
set rjL = 1 in our computation.^ By parameter rescaling, the 
above models could be connected easily to the ones gener¬ 
ally found in the literature (see e.g.. Refs. [34, 35]). 

For all three generations of LQs, we implement Model A 
in the FeynRules package [36] together with the SM in or- 


^ In general, the LHC is insensitive to this parameter. However, in case of 
a third generation LQ that couples to a f-quark, it might be accessible 
via the spin of the f-quark. 


der to collect all the tree level couplings and fields in Uni¬ 
versal FeynRules Output (UFO) format. Then, we renor¬ 
malize the Lagrangians to get rid of the UV divergences 
that appear in loop calculation. It is well known that any 
one-loop amplitude can be written as linear combinations 
of scalar tadpole, bubble, triangle and box integrals together 
with rational terms that come from the (c? — 4) part of a d- 
dimensional integral. Within the framework of the Ossola, 
Papadopoulos, and Pittau (OPP) reduction technique [37], 
the {d — 4) part that originates form the denominator of the 
integrand is called the Ri term, whereas the other rational 
term i ?2 comes due to the (d — 4) component of the nu¬ 
merator of such an integrand. We make use of the NLOCT 
package [38] to calculate the UV counterterms within the 
on-shell renormalization scheme and to determine the ra¬ 
tional i ?2 terms required for the OPP reduction. The Ri 
terms can be calculated as a four-dimensional integration 
using a particular set of scalar integrals [39]. Like Ref. 
[28], we also neglect the lepton exchange diagram (see Fig. 
1) that depends on the LQ-lepton-quark coupling (A) and 
contributes at 0{X'^) in the cross section. For the first gen¬ 
eration LQs, we estimated its contribution to be < 10% of 
the LO QCD mediated contribution for A as large as 0.5 
at the 8 TeV LHC [29]. For higher generations, it would 
be even smaller because of the relative suppression of the 
initial PDFs. Hence, all our cross section estimations are 
essentially model independent and are equally applicable 
for Model B too. An analysis for the ^ijj channel might 
also be used for a LQ with charge 5/3 (—4/3) that couples 
to a M-type (d-type) quark and a charged lepton (but not to a 
neutrino) [29] as long as £^’s are not distinguished. How¬ 
ever, only for the third generation, we implement Model 
B separately since a f-quark decays inside the detector and 
leaves very different signatures from the 6-quark. 

As already mentioned, the complete computational setup 
for the NLOh-PS correction of the LQ pair production is 
prepared under the automated MadGraph5_aMC @NLO 
environment [33]. We use MadGraph5 [40] for the 
Born level computations. The real emission corrections are 
computed with MadFKS [41] which uses the subtraction 
scheme of Frixione, Kunszt, and Signer (FKS) [42]. To 
evaluate the virtual contribution, we use MadLoop [43] 
which relies on the OPP reduction scheme. After the gen¬ 
eration of NLO events, the decay of each LQ into a lepton- 
jet pair is organized by MadSpin [44], which takes care of 
decays of heavy resonances preserving all spin information 
at the tree level accuracy. The NLO events, thus decayed, 
are then matched to the Pythia8 [45] PS following the 
MC@NLO formalism [32]. The interfacing between dif¬ 
ferent packages is made automated in aMC @NLO [46]. 

We fix the electroweak input parameters (1) Mz = 
91.188 GeV, {ii) Gf = 16.637 x 10"® GeV^ and {Hi) 
~ 132.237 in our computations. The electroweak 
mixing angle and the mass of the FU^-boson are calcu¬ 
lated from these three independent input parameters. For 
the PDFs we use MSTW(n)lo2008cl68 sets [47] that de¬ 
termine the value of strong coupling as{Mz) at (N)LO in 
QCD and use n/ = 5 massless quark flavors, i.e., the bot¬ 
tom quark is also treated effectively as a massless quark, 
since its running mass decreases significantly as the hard 
scale increases. Unless otherwise stated, we set = 650 
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MtJTeY) 


unlo (fb) 


T^-factor 


MSTW(n)lo2008cl68 PDFs 


0.2 

(7.757 

± 

0.082) 

X 

lO'* 

1.363 

0.4 

(2.295 

± 

0.019) 

X 

10® 

1.375 

0.6 

(2.292 

± 

0.018) 

X 

10® 

1.372 

0.8 

(3.829 

± 

0.033) 

X 

10® 

1.376 

1.0 

(8.190 

± 

0.065) 

X 

10° 

1.336 

1.2 

(2.113 

± 

0.015) 

X 

10° 

1.315 

1.4 

(6.017 

± 

0.046) 

X 

10-® 

1.280 

1.6 

(1.907 

± 

0.018) 

X 

10-® 

1.295 

1.8 

(6.253 

± 

0.052) 

X 

10"® 

1.292 

2.0 

(2.116 

± 

0.014) 

X 

10"® 

1.275 

— 

CTEQ6(M/L1) PDFs 

0.2 

(7.518 

± 

0.083) 

X 

10“ 

1.491 

0.4 

(2.217 

± 

0.018) 

X 

10® 

1.577 

0.6 

(2.233 

± 

0.018) 

X 

10® 

1.649 

0.8 

(3.781 

± 

0.032) 

X 

10® 

1.724 

1.0 

(8.245 

± 

0.071) 

X 

10° 

1.742 

1.2 

(2.210 

± 

0.015) 

X 

10° 

1.810 

1.4 

(6.548 

± 

0.059) 

X 

10"® 

1.852 

1.6 

(2.108 

± 

0.020) 

X 

10"® 

1.919 

1.8 

(7.113 

± 

0.053) 

X 

10"® 

1.988 

2.0 

(2.486 

± 

0.018) 

X 

10"® 

2.038 


Vs 

QCD order C -2 

C-i 

Co 

Cl 

C2 

8 TeV 

LO 

-0.209 

2.655 

2.597 

-5.534 

-0.831 


NLO 

-0.227 

2.814 

2.524 

-5.226 

-0.946 

13 TeV 

LO 

-0.283 

3.184 

2.708 

-3.931 

-0.256 


NLO 

-0.300 

3.318 

2.762 

-3.780 

-0.299 

33 TeV 

LO 

-0.217 

2.581 

6.059 

-3.617 

0.203 


NLO 

-0.221 

2.590 

6.377 

-3.621 

0.199 

100 TeV 

LO 

-0.196 

2.343 

8.613 

-3.014 

0.232 


NLO 

-0.196 

2.307 

8.909 

-2.981 

0.225 


TABLE II: Coefficients of the cross section fitting functions de¬ 
fined in Eq. (6) for different center-of-mass energies (y/S) at the 
LHC computed with MSTW(n)lo2008cl68 PDF sets [47], 


are also shown in the lower panels. The cross sections (LO 
and NLO) can be fitted with the following fitting function 
to a very good accuracy, 


a (in fb) 




exp 
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TABLE I: NLO cross sections and Ff-factors computed in two 
different PDFs for the 14 TeV LHC. 


GeV with A = 0.3. 

The central values of the renormalization and fac¬ 
torization (p_f) scales are set equal to the standard 
choice [12, 28]. To estimate the scale uncertainties of vari¬ 
ous observables, we set 

= MF = C ; (5) 

with C = {1/2,1,2}. 

Events are generated with extremely loose cuts to get un¬ 
biased results. At the time of showering these events with 
PythiaS, jets are clustered using the anti-fcx algorithm 
[48] with radius parameter r = 0.5. Before we present 
the NLOh-PS results in the next section, we note two points 
in favor of the validity of our setup: (i) we find exact nu¬ 
merical cancellation between the double and single poles 
coming from the real and virtual parts and ( ii) our FO NLO 
results are in extremely good agreement with the NLO re¬ 
sults already available [28]. 


III. NUMERICAL RESULTS 

In Table I, we compare the FO NLO cross sections 
and the iC-factors for 14 TeV LHC computed with 
CTEQ6(M/L1) PDFs [49] and n/ = 5 massless quark fla¬ 
vors (i.e., the same scheme as in Ref. [28]) with those ob¬ 
tained with our scheme. The numbers obtained with CTEQ 
PDFs can be compared directly with the ones shown in Ta¬ 
ble I of Ref. [28]. 

In Fig. 2, we show the LQ pair production (pp iqiq) 
cross sections computed with C = 1 [Eq. (5)] for four dif¬ 
ferent center-of-mass energies {\fS) at the LHC environ¬ 
ment (8,13, 33 and 100 TeV). The corresponding AT-factors 


where the coefficients are given in Table II. These fits are 
valid for between 0.25 and 4 TeV for all the center-of- 
mass energies except y/S = 8 TeV for which the range of 
validity is 0.25 TeV < Mi^ < 2 TeV. 

To show the effect of the PS over the NLO FO 
results, we plot the transverse momentum distribution, 
LQ-pair in Fig. 3 for the 14 TeV 
LHC. The fixed order (NLO FO) distribution is shown by 
the dashed blue line, and the NLO matched with the par- 
ton shower (NLOh-PS) result is shown with a solid brown 
line. The NLO FO result is divergent as px —5" 0 and there¬ 
fore it is necessary to match it with PS in order to get a 
reliable estimation of the cross section in this region. The 
converging behavior of the NLOh-PS distribution in the low 
px region indicates the correct resummation of the Sudakov 
logarithms in the collinear region, thereby leading to a no¬ 
table Sudakov supression visible prominently in the log-log 
plot of Fig. 3. As expected, both of these results agree in 
the high px region. 

Since the analyses for the first two generations of LQs 
are very similar, in the next two figures, we present results 
of kinematic distributions of a selection of observables at 
(N)LOh-PS accuracy only for the first generation. Note that, 
from here onward, we set the typical value of LHC center- 
of-mass energy to yfS =13 TeV unless specified otherwise 
and pass the Pythia events through the detector simulator 
Delphes 3.3.1 [50] with the default CMS card to gener¬ 
ate the distributions. To obtain the distributions, we apply 
the basic preselection kinematic cuts from the CMS analy¬ 
sis^ [12] to the eejj final states: 

(i) exactly two electrons (e^) with transverse momen¬ 
tum p^ > 45 GeV and pseudorapidty \rje\ < 2.5 


^ These cuts are designed for the 8 TeV LHC. However, since they are 
only for the preselection, we use them to get the 13 TeV distributions 
for easy comparison. 
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FIG. 2: The LO and NLO cross sections computed with MSTW(n)lo2008cl68 PDFs [47] and fip ~ = Mi, for four different 

center-of-mass energies at the LHC environment - (a) 8 TeV, (b) 13 TeV, (c) 33 TeV and (d) 100 TeV. The corresponding TT-factors are 
shown in the lower panels. The lines represent the fitting functions defined in Eq. (6) with the coefficients given in Table II whereas the 
points are obtained from direct computations. 


Transverse momentum of LQ-pair 



FIG. 3: Transverse momentum distribution of the LQ-pair at fixed 
order NLO and NLO-l-PS at the 14 TeV LHC. 


excluding 1.442 < |? 7 e| < 1.56, 

(ii) two leading pT jets with > 125 GeV, > 45 
GeV and \r]j\ < 2.4, 

(Hi) separation between an electron and a jet in the t] — (f> 
plane, ARej > 0.3 

(iv) the invariant mass of the electron pair, Mge > 50 
GeV, 


(v) the scalar sum of the px of the two electrons and the 
two leading px jets, St > 300 GeV, 

(vi) the minimum of electron-jet invariant mass com¬ 
binations, M™'" > 50 GeV, where, to form two 
ej-pairs out of eejj final states, combination with 
the smaller difference between the two ej invariant 
masses (Mej) is considered. 

Similarly, for the ei/jj final states, we apply 

(i) exactly one electron with p^ > 45 GeV and |r7e| < 
2.1 excluding 1.442 < |pe| < 1.56, 

(ii) two leading px jets with p?^ > 125 GeV, pi^ > 45 
GeV and \rij\ < 2.4, 

(Hi) the missing transverse energy, > 55 GeV, 

(iv) azimuthal separation between the electron and the 

> 0 . 8 , 

(v) azimuthal separation between the hardest jet and the 

> 0.5, 

(vi) separation between the electron with the two leading 
Px jets in the r] — (p plane, AR^j > 0.7, 
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Average electron-jet invariant mass 



Minimum electron-jet invariant mass 



(e) 


(f) 


FIG. 4: Kinematic distributions for the eejj channel. For these plots we set Mi^ = 650 GeV and /3 = 1. 
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(a) 


(b) 





(e) (f) 


FIG. 5: Kinematic distributions for the evjj channel. For these plots we set = 650 GeV and /3 = 0.5. 
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FIG. 6: Distributions for third generation charge 1/3 LQs. For these plots, we set = 650 GeV and /? = 0.5. See Eq. (8) for the 
definition of Mr. 


(vii) the scalar sum of the px of the electron, the 
and the two leading jets, St > 300 GeV, 

(viii) the electron-neutrino transverse mass, = 

[1 - cos (A(/> e))] > 50 GeV, 

(ix) the electron-jet invariant mass, Mej > 50 GeV, 
in which the electron and the jet are picked out of 
that combination for which the difference between 
the electron-jet transverse mass and the neutrino-jet 
transverse mass is smaller. 

To present the plots, we adopt a systematic graphical rep¬ 
resentation scheme; the main frame shows the kinematical 
distributions of an observable at LOh-PS (dashed blue) and 
NLOh-PS (solid brown) accuracy, the middle pad represents 
the ratio (7f-factor) of the NLOh-PS result over the LOh-PS 
result with a solid black line and in the lower inset, and we 
present the fractional scale uncertainties of corresponding 
observables [dashed blue (solid brown) for (N)LOh-PS], de¬ 
fined by the ratio of their variation around the central scale. 

In Fig. 4, we show the distributions for the eejj chan¬ 
nel. Our choice of variables is mainly motivated by 
the CMS analysis [12] (ALTAS also uses similar vari¬ 
ables [6]). In addition to the px of the hardest elec¬ 
tron and jet [Figs. 4a and 4b], we show the distribu¬ 
tions of St, Mee, (the average of electron-jet in¬ 

variant mass combinations) and M™'" in Figs. 4c, 4d, 
4e and 4f respectively. For the eiyjj channel, we show 
the distributions in Fig. 5: (i) [Fig. 5a], (ii) 

[Fig. 5b], (iii) St [Fig. 5c], (iv) M,, [Fig. 5d], (v) 
MT,eii [Fig- 5e], and {vi) the neutrino-jet transverse mass 

MT,vj = \j (Lix*’*® + E^)'^ - +Px)^ [Fig. 5f] 

where the combination with the smaller difference between 
the electron-jet transverse mass and the neutrino-jet trans¬ 
verse mass was considered. 

For all the kinematical observables presented in Figs. 4 
and 5, we observe that the NLOh-PS distributions are sub¬ 
stantially larger than the corresponding LOh-PS results and 


the scale uncertainties are notably reduced at NLOh-PS. 
Though not shown, we find no significant differences be¬ 
tween the PDF uncertainties computed at the NLOh-PS and 
LOh-PS levels using Hessian method as suggested by the 
MSTW Collaboration [47]. 

In Figs. 4a and 4b or 5a and 5b we find that in the high 
px/ii/xregions, there is no notable difference between 
NLOh-PS and LOh-PS; i.e., the AT-factor is close to 1 in 
those regions. This is because, while the or the px of 
an electron or the leading jet increases, the extra emission at 
NLO tends to have smaller px as the center of mass energy 
remains unaltered. A similar argument holds for Fig. 4c 
or. 5c. In the other distributions, the AT-factor is nearly 
constant in the regions shown; however with the increase 
in mass scale it comes down to 1 as it should. 

For the third generation, as already mentioned, we have 
implemented both Model A and Model B in our code. Here, 
however, we only display the plots for Model A type LQs 
that decay either to f-quarks and r’s or to 6-quarks and i/^-’s. 
In Fig. 6a, we show the number of hard (with px > 50 
GeV) electrons or muons in the events generated with /3 = 
0.5. Notice, in this case there can be more than two charged 
leptons in the final states from decays like 

r—decay 

fr (6 f {i Vt (7) 

decay 

where the two light leptons (e or p) coming from the same 
LQ will have opposite charges. For a LQ with electric 
charge ±4/3, they would carry the same charge but such a 
LQ would not couple to a z/ and a 6-quark i.e., /3 cannot be 
0.5 (as we have set here), in the absence of any unknown 
decay mode. This multilepton signature is unique for the 
third generation, as pair productions of the first two genera¬ 
tions of LQs can produce only two hard leptons in the final 
state. 

For the events where the LQs have decayed into 6-quarks 




















8 


and Vr’s we have plotted Mu (see Fig. 6b), defined as 

Mr = (l^p + 1^1)^ - (pi (8) 

where p* refers to the ith px ordered jet and p* its lon¬ 
gitudinal component. It is an approximation of the razor 
mass used in Ref. [11], as we have used the first two hard 
jets instead of “pseudojets” (see Ref. [11] for details). It is 
suitable for identifying a heavy resonance decaying into a 
jet and a weakly interacting neutral particle from the SM 
background. For our demonstration plot, we have selected 
events with no light charged lepton and >140 GeV. 

IV. CONCLUSIONS 

LQs were hypothesized a few decades back but mainly 
because of the failures of earlier LQ-searches at HERA and 
Tevatron, nowadays they receive much less attention than 
other BSM candidates. However, they are still relevant to¬ 
day (see e.g.. Refs. [51-53]) and the LHC is actively look¬ 
ing for their signatures. 

In this paper, we have taken another look at the NLO 
QCD corrections to the scalar LQ pair production process 
at the LHC and included the parton shower effects to obtain 
new NLOh-PS accurate results. Our computation provides 
a reliable and completely automated code that goes beyond 
computing the NLO iC-factor by simulating signal events 
at the QCD NLO(h-PS) level and improves the precision 
of various kinematic distributions. These precise distribu¬ 
tions can play a crucial role in advanced techniques like 
multivariate analysis. The excellent agreement between 
our cross section computations with the available results 


demonstrates the reliability of our code. It is prepared us¬ 
ing the MadGraph5_aMC@NLO framework and there¬ 
fore, is flexible to the choice of energy, cuts, scales, PDFs, 
etc. We stress that our computations cover all three genera¬ 
tions. The signatures of the pair productions of the first two 
generations of LQs are similar and mostly model indepen¬ 
dent, whereas for the third generation, they are different in 
nature. Our code can be used for LQs decaying to either a 
top quark or a bottom quark in association with a third gen¬ 
eration lepton. The complete stand-alone code is available 
publicly at the website - http: //amcatnlo. cern. ch. In 
addition, we have also provided A-factors for a wide range 
of Mf,^ by means of fitting functions for four different LHC 
center-of-mass energies, ^fS =8, 13, 33, 100 TeV. 

Finally, we note that sometimes searches for other BSM 
particles with similar final states are reinterpreted for LQ 
parameters. For example, in Ref. [ 11 ], a search for particles 
decaying to a 6-quark + A™'®® final state is interpreted in 
terms of 6-squarks decaying to 6-quarks and xo’s and third 
generation LQs decaying to 6-quarks and z/t-’s. However, 
this could be misleading if the xo has large mass. With our 
code, such problems can be avoided by simulating the LQ 
signal precisely and independently. 
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